Example Usage

%matplotlib inline
%config InlineBackend.figure_formats = ['svg']
import chalc as ch
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["animation.html"] = "jshtml"
import h5py, warnings, logging, sys

# Disable some deprecation warnings from the seaborn plotting library
logging.basicConfig(stream=sys.stdout, level=logging.ERROR)
logging.captureWarnings(True)

For our data we sample 100 points on a circle with some noise and 100 points from inside the unit disk.

rng = np.random.default_rng(40)
num_points_circle = 100
num_points_disk = 100
mean = [0, 0]
cov = np.eye(2)*0.01
circle = np.array([[np.sin(2*np.pi*t), np.cos(2*np.pi*t)] for t in rng.random(num_points_circle)]).T +\
    rng.multivariate_normal(mean, cov, num_points_circle).T # points as columns
disk = np.sqrt(rng.random(num_points_disk)) * np.array([[np.sin(2*np.pi*t), np.cos(2*np.pi*t)] for t in rng.random(num_points_disk)]).T
points = np.concatenate((circle, disk), axis=1)
plt.scatter(circle[0, :], circle[1, :], s=10)
plt.scatter(disk[0, :], disk[1, :], s=10)
plt.gca().set_aspect('equal')
colours = [0]*num_points_circle + [1]*num_points_disk

We compute the chromatic alpha complex \(K\) of the point cloud. \(K\) has far fewer simplices than either the Cech or Vietoris-Rips complex, which each have \(\displaystyle \binom{200}{2} = 19900\) edges and \(\displaystyle \binom{200}{3} = 1313400\) 2-simplices.

K = ch.chromatic.alpha(points, colours)
print(f'{len(K.simplices[1])} 1-simplices')
print(f'{len(K.simplices[2])} 2-simplices')

954 1-simplices
1312 2-simplices

We can visualise the filtration at any given time. For example, at time 0.2 we have:

ch.plotting.draw_filtration(K, points, 0.2)

(<Figure size 640x480 with 1 Axes>, <Axes: >)

We can compute the six-pack of persistent homology diagrams of the chromatic alpha complex, the chromatic Delaunay-Cech complex, and the chromatic Delaunay-Rips complex.

dgms_alpha = ch.sixpack.compute(points, colours, dom=[0,], method="chromatic alpha")
dgms_delcech = ch.sixpack.compute(points, colours, dom=[0,], method="chromatic delcech")
dgms_delrips = ch.sixpack.compute(points, colours, dom=[0,], method="chromatic delrips")

Each six-pack is an object having attributes ker, cok, im, dom, cod, rel, entrance_times, and dimensions. Each of the first six of those attributes is an object that has sets of paired and unpaired simplices, while the last two attributes indicate the filtration time and dimension of those simplices. For example-

print(dgms_alpha.ker.paired)

{(1187, 1200), (1445, 1453), (358, 409), (1995, 1999), (584, 883), (1697, 1701), (228, 235), (333, 581), (361, 456), (454, 471), (1648, 1650), (1486, 1499), (1447, 1495), (1461, 1494), (1565, 1600), (429, 502), (1510, 1514), (606, 626), (1360, 1458), (1248, 1252), (1393, 1403), (1267, 1268), (487, 531), (387, 587), (506, 596), (1454, 1457), (292, 383), (348, 528), (1329, 1331), (421, 491), (252, 438), (1592, 1671), (280, 317), (1464, 1470), (1366, 1408), (359, 412), (2013, 2433), (345, 376), (282, 338), (413, 428), (544, 646), (352, 540), (262, 372), (1692, 1693), (436, 658), (1505, 1516), (331, 416)}

A convenience function is provided to extract the diagrams as a matrix of birth/death times:

np.set_printoptions(threshold=10)
print(dgms_alpha.get('ker', 1)) # get the kernel in dimension one
# dgms_alpha.get('ker', [0, 1]) to get the kernel in dimensions zero and one
# dgms_alpha.get('ker') to get the kernel in all dimensions from zero to max(dgms_alpha.dimensions)

[[0.0566658  0.06077279]
 [0.10478006 0.10515716]
 [0.22895679 0.23057114]
 ...
 [0.23772323 0.78357409]
 [0.14216739 0.14226804]
 [0.1135806  0.1155313 ]]

Convenience functions are provided for plotting either a six-pack or a specific diagram from a six-pack.

fig1, ax1 = ch.plotting.plot_sixpack(dgms_alpha)
fig1.suptitle('Chromatic Alpha')

fig2, ax2 = ch.plotting.plot_diagram(dgms_delcech, 'rel')
fig2.suptitle('Chromatic Delaunay-Cech')

Text(0.5, 0.98, 'Chromatic Delaunay-Cech')

We can save the diagrams to file or load a diagram from file:

with h5py.File('test.h5', 'w') as f:
    ch.sixpack.save_diagrams(dgms_alpha, f)

with h5py.File('test.h5', 'r') as f:
    dgms_alpha = ch.sixpack.load_diagrams(f)

We can visualise the 2-skeleton of the filtration for points in 2D:

animation = ch.plotting.animate_filtration(
    K, points, filtration_times=np.linspace(0, 1.0, 45), animation_length=5)
animation

Once Loop Reflect